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We present a study of the solvation properties of model aqueous electrode interfaces. The exposed 
electrodes we study strongly bind water and have closed packed crystalline surfaces, which template 
an ordered water adlayer adjacent to the interface. We find that these ordered water structures fa- 
cilitate collective responses in the presence of solutes that are correlated over large lengthscales and 
across long timescales. Specifically, we show that the liquid water adjacent to the ordered adlayers 
forms a soft, liquid- vapor-like interface with concomitant manifestations of hydrophobicity. Tem- 
poral defects in the adlayer configurations create a dynamic heterogeneity in the degree to which 
different regions of the interface attract hydrophobic species. The structure and heterogeneous dy- 
namics of the adlayer defects depend upon the geometry of the underlying ordered metal surface. 
For both 100 and 111 surfaces, the dynamical heterogeneity relaxes on times longer than nanosec- 
onds. Along with analyzing time scales associated with these effects, we highlight implications for 
electrolysis and the particular catalytic efficiency of platinum. 



Extended metal interfaces play a fundamental role in 
aqueous electrochemistry, a field of principle importance 
in the advancement of renewable, clean energy sources 
Pfl [2]. In many processes that occur at metal inter- 
faces, such as electrolysis, corrosion and electrocataly- 
sis, water is ubiquitous, often acting as both solvent and 
reactant[3]. While many studies exist detailing the be- 
havior of water across small length and timescales, [1H7] 
and at low temperatures [HI [9], at present there is lit- 
tle understanding of the large lengthscale correlations 
and emergent behavior of water on metal surfaces, even 
though such effects are likely to influence function in im- 
portant ways [T0UT3] . Here, we report classical simula- 
tion results for water on an atomistic model electrode 
where a detailed study of solvation at the surface has 
been undertaken. In this study we place particular em- 
phasis on long time correlations and large lengthscale 
phenomena. Specifically we illustrate how an electrode 
can impose geometrical constraints within the adlayer of 
water, creating a composite metal-water interface that 
is hydrophobic on large length scales. We further show 
how defects within the hydrogen bonding patterns of the 
adlayer create transient regions of hydrophobic behav- 
ior that exist on small length scales and over long times 
scales. These results offer a microscopic explanation and 
generalization of previous experimental observations that 
have inferred surface hydrophobicity of a platinum elec- 
trode at low temperatures [14] [15] . 

To study the electrode interface we use molecular mod- 
els [6] [16] that neglect explicit electronic degrees of free- 
dom beyond accounting for electronic polarization of the 
electrode. Despite its relative simplicity, the model shows 
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reasonable agreement with experiment for the potential 
of zero charge and capacitance [10 of the aqueous plat- 
inum interface. For example, we can estimate the po- 
tential of zero charge, U pzc relative to the hydrogen elec- 
trode following previous work [17] by measuring the con- 
tact potential between the electrode and the water bi- 
layer, ip w —0.9V, and referencing it to the known work 
function for platinum, W = 5.9V [17] and the absolute 
hydrogen electrode potential, Un 2 = 4.4V. Using this ap- 
proach we get U pzc = 0.6V compared to the experimental 




FIG. 1. Illustration of the separation of timescales between 
reorganizing surface configurations, which occur on average 
every r s , and reorganizing the bulk density, which occur on 
average every Tt>. Thin tic marks are separated by 20 ps, 
which is on the order of though larger than timescales for 
typical density fluctuations. Thick tick marks are separated 
by 100 ps, which is of the order of the typical relaxation times 
for relevant interfacial fluctuations. 
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value of U pzc = 0.4V[T8]. 

Using this model we can access length and time scales 
far beyond those currently available to ab inito calcula- 
tions, and probe collective behavior. We find that the hy- 
drophobicity of the water-electrode interface depends on 
the amount of passivation of the hydrogen bond network 
within the adsorbed water layer. Strong, favorable inter- 
actions of the order of half an eV between the water and 
electrode pin the oxygens of the water to the top sites of 
the crystal lattice creating a spatially ordered arrange- 
ment of molecules. We have studied two surfaces of a 
planar electrode, corresponding to the 100 and 111 faces 
of an FCC crystal, with lattice spacings corresponding 
to a platinum lattice. On both the 100 and 111 surfaces, 
the imposed water structures allow for facile hydrogen 
bonding within the adlayer and subsequently only a few, 
fleeting, hydrogen bonds are donated from the adlayer 
to the surrounding bulk. While the adsorbed oxygens at 
both surfaces still afford hydrogen bond acceptor sites, 
the asymmetry associated with lacking donor sites re- 
sults in an interface that is liquid- vapor like in the sense 
that large density fluctuations occur though the collective 
formation and deformation of an interface [15]. 

Even though the underlying electrode lattices we study 
are perfectly ordered, over large length scales the pla- 
nar geometry of the surface is incommensurate with wa- 
ter's preferred tetrahedral structure. A consequence of 
this frustration is the presence of an equilibrium num- 
ber of defects in the hydrogen bond network within the 
adlayer. These defects facilitate reorganization within 
the surface, as detailed elsewhere [20]. and the resulting 
dynamics are heterogeneous and relax on timescales of 
1-10's of nanoseconds. This timescale for surface relax- 
ation, which we denote r s , is much larger then typical 
times for equilibrium density fluctuations, r s ^> Tb « 5ps. 
Therefore, while the presence of the surface introduces a 
static inhomogeneity into the system, the water bound to 
this surface introduces a dynamic inhomogeneity into the 
system. The resultant separation of timescales between 
bulk and surface reorganization is shown in Fig. [I] The 
top panels of Fig|l] show snapshots during slow reorga- 
nization of the surface water dipoles, while the bottom 
panels of Fig [I] illustrate faster interfacial fluctuations. 

These two features of the electrode interface, the static 
heterogeneity of the extended interface and the slow dy- 
namics of water at its surface, cause a decoupling of en- 
semble and dynamic averaging on timescales t < r s . For 
a given configuration of the adlayer the liquid water in 
contact with it swiftly equilibrates and for t ^> the 
properties of the subsequent solvation layers are station- 
ary. Over timescales intermediate between surface and 
bulk relaxation, we find that the temporal heterogeneity 
of the hydrogen bond network couples to the solvation 
properties of the interface, creating transient regions of 
favorable and unfavorable solvation. These regions of dif- 
fering solvation interconvert on the timescales of surface 
relaxation, and only for t > r s , does the solvation at the 
interface reflect the homogenous symmetry of the fluid 



above it. 

The rest of the paper expands on these findings. First, 
we illustrate the particular structure of the adlayer and 
show how the passivation of the hydrogen bond network 
on the surface creates a liquid-vapor-like interface that 
easily accommodates small, spherical and large, extended 
hydrophobes. Next, we show how the strain in the wa- 
ter structure on the surface, coupled with a separation 
of relaxation times between the surface and bulk, cre- 
ates temporal regions of spatially heterogeneous solva- 
tion potential that decay over nanoseconds. Finally we 
discuss how the presence of a liquid-vapor like interface 
is expected to influence electrochemical properties by at- 
tracting specific ions such as excess protons necessary for 
chemical reactions. 



I. WATER ELECTRODE INTERFACE 

The properties of the electrochemical interface are dic- 
tated in large part by the first adsorbed water layer. The 
generic features of this surface are that 1) single water 
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FIG. 2. Illustration of the 100 and 111 adlayers and their 
effect on macroscopic solvation, (a) The 100 surface is locally 
four coordinated and commensurate with favorable hydrogen 
bonding patterns. Large ordered domains are separated by 
line defects, (b) The highly ordered domains donate few hy- 
drogen bonds to the subsequent water layers leading those 
layers to not wet the composite water electrode surface, (c) 
The 111 surface is locally six coordinated and frustrates hy- 
drogen bonding. Though water is still ordered, vacancies and 
interstitial defects are common, (d) Hydrogen bond donor 
sites are more common strongly than on the 100 surface, and 
subsequently the surface is wet. 
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FIG. 3. Structure and solvation of the composite water-electrode interface, (a) The density profile of water molecules away 
from both surfaces, divided by the bulk water density, p w — 0.033A -3 . The density is characterized by a strongly bound 
adsorbed layer and solvent layering that extends roughly 10 A into the bulk. Immediately adjacent to the adsorbed layer is a 
region of density depletion, (b) For the 100 surface, shown, or a 111 surface, a typical configuration of water molecules shows 
very little hydrogen bonding (dashed red lines) between the surface and the bulk. The volume v (black circle) is a 3 A radius 
sphere, (c) The excess solvation free energy for a 3 A in radius ideal hydrophobe, shown to scale in (b), as a function of 
distance away from the surfaces. The dashed lines in (a) and (c) define z* as the distance of closest approach of the 3 A sphere 
to the adlayer. 



molecules bind strongly and specifically to the top sites 
of the metal and 2) the crystalline geometries of these top 
sites are incommensurate with extended hydrogen bond- 
ing patterns. Quantum chemical calculations [5] and ex- 
perimental photoelectron spectroscopy measurements [9 
have estimated the single molecule binding energy of wa- 
ter a platinum electrode to be close to 0.4 eV, with a 
small dependence on surface geometry. The model used 
here [6 has been parameterized to recover this strong 
attraction with binding energies of 0.46 eV and 0.37 eV 
for the 100 and 111 surfaces respectively. Depending on 
the specific surface geometry, the thermodynamically sta- 
ble structure may not be a complete monolayer. This is 
due to the competition between surface adsorption and 
hydrogen bond formation within the surface, which is 
frustrated on surfaces incommensurate with water's pre- 
ferred hydrogen bonding network which favors four fold 
coordination [2T]. 

This interplay of energy scales is reflected in the struc- 
tural motifs present on the different crystal faces. Fig- 
ure [2] shows characteristic snapshots of the adlayer of wa- 
ter for both surfaces (a and c) as well as their subsequent 
effect on wetting (b and d). For the 100 surface, metal 
atoms are locally four fold coordinated and are thus mod- 
erately commensurate with a two-dimensional projection 
of local hydrogen bonding patterns. As a result, the 
structure of water on the surface is highly ordered with 
the water's dipole oriented parallel to the surface and ap- 
proximately all top sites are occupied. At any particular 
instant, however, line defects exist on the surface sep- 
arating planes of dipole aligned molecules by 90 degree 
turns in their orientations. For the 111 surface, metal 
atoms are locally six fold coordinated and though they 
also have lattice spacings that are commensurate with a 



hydrogen bond, the six-fold coordination frustrates pre- 
ferred bonding patterns. As a result this surface has 
regions of local hexagonal order, rings are seen in the 
structures adopted by monolayers of water absorbed on 
the 111 surface of many face-centered cubic metals such 
as Pt and Pd [22 , together with an equilibrium concen- 
tration of interstitials that occupy empty top sites with 
dipoles that point away from the surface. This disorder 
results in an average coverage of about 85% of all top 
sites. For both surfaces the lattices are entirely regular, 
and therefore the heterogeneity in the hydrogen bond- 
ing network is temporary. However the imposed order 
within the adlayer dictates that relaxation occurs over 
long timescales[20]. 

The presence of extended interfaces in solution, such as 
the solvated electrode surface, are expected to influence 
the properties of subsequent solvent layers over distances 
corresponding to the bulk correlation length. For a liquid 
near coexistence with its vapor, such as water at ambient 
conditions, extended inhomogeneities can give rise to a 
de- wetting transition [19], whose interfaces subsequently 
have larger correlation lengths. Figure [3ja) illustrates 
the mean density of water molecules as a functions of 
the distance away from the electrode surface. Although 
the structure on the adlayer depends intimately on the 
electrode geometry, the surrounding water is fairly in- 
sensitive to the exposed crystal face. We find for both 
surface geometries that the density profile for water away 
from the interface exhibits a sharp peak at the electrode 
surface, indicative of the adlayer, followed by a region 
of a density depletion approximately 3 A thick. Density 
oscillations decay over 1 nm away from the surface. The 
region of depleted density demonstrates that the asym- 
metry between hydrogen bond donors and acceptors at 
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the interface results in an unbalanced attraction, how- 
ever the effective interaction with the surface is not so 
weak so as to allow the formation of capillary waves that 
would destroy the density oscillations seen away from the 
electrode. 

The density depletion seen immediately adjacent to the 
adlayer is enough to make solvation of ideal hydrophobes 
(hard-spheres) favorable at the interface. Figure [3] (b) il- 
lustrates the size of a 3 A in radius sphere positioned near 
the interface. In this representative snapshot it is clear 
that typical hydrogen bonding patterns near the interface 
easily accommodate the solvation of such probe volumes. 
Figure [3] (c), illustrates the excess chemical potential for 
a hard sphere with a radius of 3 A. We calculate this 
quantity by monitoring the number fluctuations within a 
probe volume. Specifically, we calculate the probability 
of observing N molecules in the probe volume, u, 

P V (N) = (5(N V - N)} 
1 ( l 

= lim - / dt'5(N v (t') - AT), (1) 
t^oo t J 

where S(N V — N) is a Kroniker delta function and (...) 
denotes equilibrium average. This distribution is re- 
lated to the excess solvation free energy for an ideal hy- 
drophobe through the relation [23 



Pv(0) 



(2) 



where Aja v is the reversible work to create a cavity of 
size and shape v and f3 is Boltzman's constant. 

For the systems we consider here, the existence of the 
planar electrode breaks transitional invariance. In order 
to accommodate this aspect we denote, P V ^(N), where 
r is the position of the center of the probe volume. This 
distribution reduces to P v ( r \(N) — » P V (N) when r is far 
away from the electrode surface. Correspondingly, we 
also define P v ^(0) = exp [— /3A/i v ( r )]. In other words, 
solvation free energy in an inhomogeneous system is gen- 
erally spatially dependent. 

Due to the separation of timescales between surface 
and bulk relaxation, our system is also dynamically 
heterogeneous. Therefore, on intermediate timescales, 
Tb << t < r s , the solvation free energy carries a time- 
dependence. This time dependence is denoted as, 



1 



dt f 5[N v(r) (t f ;x )-N] (3) 



_ e -/3A/x„( r )(t; x ) 



(4) 



where xo denotes the initial surface configuration and t is 
the timescale over which the distribution is averaged. For 
t <C t s , Eq. E simplifies to, P v ^(N,t; xo) — » P V ^(N). 
Finally, the difference between the value of the solvation 
free energy located at r, averaged over a time t, and its 
equilibrium bulk value is defined as 



For t — )► 00, S/jl v (t, t; x ) — » Sfi v (r) and for r far away 
from the surface 5[i v (r) —> 0. At long times, and av- 
eraged over the plane of the surface, the solvation free 
energy has a minimum at the distance of closest ap- 
proach to the composite water-electrode surface, z*, in- 
dicated by a dashed line in Fig. [3J The negative sol- 
vation free energy implies that while the bare electrode 
surface attracts water, the composite water-electrode sur- 
face preferentially attracts oil and can be considered hy- 
drophobic. While both electrode geometries exhibit en- 
hanced hydrophobic solubility, the 100 surface is more 
hydrophobic as measured by its excess solvation free en- 




FIG. 4. Solvation of a large probe volume near the interface, 
(a) A configuration of water molecules near a 100 surface with 
the instantaneous liquid interface shown in blue 24]. This con- 
figuration is typical of one solvating a cuboid of size 20 by 20 
by 3 A, shown to scale in black, (b) Probability distribution 
for finding N particles in a probe volume whose outer edge is 
located at z* t with a mean occupancy N. N = 31,34 and 40 
for the 100 surface, 111 surface and bulk, respectively. 



£/x v (r,t; x ) = A/i v(r )(t; x ) - Afi v . 



(5) 
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ergy at z* , f35n(z*) « —2.0, compared to the 111 surface, 
f35n(z*) « -0.7 for the 3 A sphere. 



II. LARGE LENGTHSCALE SOLVATION 

Using the method of indirect umbrella sampling (IN- 
DUS) [25 we are able to compute stationary distribution 
functions, P v ( r )(iV), for extremely rare fluctuations. By 
studying the tails of these distributions we can determine 
to what extent interface formation, as opposed to Gaus- 
sian density fluctuations, are important. 

The specific dimensions of the probe volume are cho- 
sen to focus on interfacial fluctuations. In particular, the 
probe volume is thin enough, 3 A, to include molecules 
that can be part of a liquid interface while not also con- 
taining molecules that are part of the bulk; and it is wide 
enougn, 20 x 20A" , to capture nano-scale fluctuations 
intrinsic to the soft liquid interface. Figure (a) illus- 
trates a configuration of water and the instantaneous liq- 
uid interface [24] solvating the cuboid sub- volume, high- 
lighting how solvation at the surface occurs by deforming 
a soft interface. The signature of this mode is illustrated 
in Figj4] (b) where highly correlated behavior, interface 
formation in this case, is apparent by the appearance of 
a exponential tail for small N for both surfaces relative 
to the bulk. These distributions along with Eqs. 2 and 
4 allow us to calculate the excess solvation free energy. 
We find that this is negative at both surfaces, however 
solvation for this large probe volume at the 100 surface 
is more favorable by 40/cbT compared to the 111 sur- 
face. While both surfaces afford hydrogen bond accep- 
tor sites, only the 111 surface has a nonzero number of 
fleeting hydrogen bond donors. Using a standard criteria 
for defining a hydrogen bond [26]. we calculate an aver- 
age hydrogen bond donor density on the surface to be 
approximately 1.0/nm 2 for the 111 surface and 0.0 /nm 2 
for the 100 surface. This means that in the large probe 
volume, there is on average 4 hydrogen bonds donated to 
the bulk. These few hydrogen bonds, produce the 40 k^T 
change in solvation, and it is expected that the addition 
of further hydrogen bond defects would make the surface 
hydrophilic. 

This microscopic measure of solvation explains the dif- 
ferent wetting behaviors observed for these surfaces in 
Fig. [2] Large lengthscale solvation at hydrophobic sur- 
faces is dominated by deforming existing interfaces, the 
excess chemical potential is expected to be well approxi- 
mated by 

/^/i(z*) = -A 7LV (l-cos#) (6) 

where A is the cross-sectional area of a large probe 
volume, 7lv is the liquid-vapor surface tension, and is 
the water droplet contact angle on the surface. The more 
favorable solvation at the surface of the 100 surface is ex- 
pected to result in a contact angle « 90°, whereas the 
subtly favorable solvation at the 111 surface is expected 
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FIG. 5. Heterogeneous solvation at the electrode surface. 
(a,b) The excess free energy for a 3 A ideal hydrophobe is 
heterogeneous at both 111 and 100 surfaces. Regions of fa- 
vorable and unfavorble solvation have been determined by 
averaging over 1 ns from an initial surface configuration, xq. 



to result in a contact angle of « 40°. These are in quali- 
tative agreement with the configurations shown in Fig. [2] 
The connection with macroscopic wetting behavior pre- 
sented here is also in agreement of previous observations 
of hydrophobicity of the adlayer for the platinum 111 
surface [14], and suggests that studies probing the plat- 
inum 100 surface would find an even more pronounced 
effect. 
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t = 1 ns ^ t = 2 us c t = 9.5 ns 




FIG. 6. Evolution of the solvation free energy as the surface reorganizes from an initial configuration, xo. (a-c) Maps of the 
excess free energy averaged over different observation times for the 111 surface. Contour plots are reported using the same 
scale shown in the colorbar on the left, (d) The time correlation function for solvation heterogeneities for each surface. Points 
(a), (b) and (c) coincides with times for panels (a), (b) and (c). 



III. HETEROGENEOUS SURFACE SOLVATION 

The data presented in Fig. [3] illustrates equilibrium 
values of the excess solvation free energy averaged over 
the plane parallel to the surface, computed by averag- 
ing over long molecular dynamics trajectories, of roughly 
10 ns. The presence of a region of strong water density 
depletion induced by the local structure of the water ad- 
layer implies that only on timescales much longer than 
the correlation time for typical bulk density fluctuations, 
Tb « 5ps [27], will solvation within this plane should be 
homogeneous, i.e. /35/j,(x, y) ~ const. However, the or- 
dering within this surface adlayer makes reorganization 
difficult, and as a result the timescale associated with de- 
correlating surface configurations, r s w 1 — 10ns, is long 
(see Appendix A). For intermediate times, Tb <C t < r s , 
this long time surface relaxation couples to the solvation 
in interesting ways. In particular, we find that for av- 
erages computed over this time, 1 ns < t < 10 ns, the 
solvation calculated within the surface is heterogeneous. 

Figure [5] depicts the solvation free energy, 
PSfiUx, y, z*}, t; xo) at both surface geometries for 
a 3 A sphere, as a function of position in the x, y 
plane and observation time, for a given surface config- 
uration. This method of spatially resolving the local 
hydrophobicity is a time dependent extension of previous 
work by others on protein surfaces, which have static 
heterogeneity [28] . As shown in Fig. [5] the structure of the 
solvation on this surface reflects neither the underlying 
lattice symmetry nor the homogeneous symmetry of the 
above liquid, but rather the coupling between hydrogen 
bonding defect structures of the bound water adlayer 
and the above solvent layers. We can quantify the 
surface heterogeneity by calculating a time-dependent 
variance 

C(t) = {[5fi({x, y, z*}, t; x ) - Sfi(z*)} 2 ) (7) 

where 5fi(z*) is an average excess solvation free energy at 
and as before (...) denote averages over realizations 
of initial surface conditions. We find for times satisfying 



t ^> Tb, the spatial average over the surface is equal to 
the long time average. Using this measure we find that 
for all times the solvation on the 111 is more heteroge- 
neous than on the 100 surface, owing to the larger domain 
sizes seen on the ordered bound layer in the 100 surface. 
These domains sizes are on average approximately 1A 
for the 111 surface and 3 A for the 100 surface, as ob- 
tained by coarse-graining Fig. [5|a,b) over 1 k^T. Figure 
[6] (a-c) show (35fi({x, y, z*}, t; xo) as t is increased for the 
111 surface. Generically, reorganization on the surface 
occurs as t is increased, and the amount of heterogeneity 
is reduced. Similar behavior is found for the 100 sur- 
face. In order to quantify the timescales for relaxing this 
heterogeneity we measure the decay of C(t)/C(0), where 
the argument of the denominator is taken at the smallest 
t where r b <C t « 0.5ns. Figure [6] (d) plots C(t)/C(0) 
for both surfaces and illustrates that the time to reach a 
uniform solvation potential at the surface is on the order 
of 10 ns. This time is on the order of many of the slow 
processes the occur on the electrode surface such as the 
mean dipole correlation time 1-10 ns (see Appendix A). 

IV. IMPLICATIONS FOR ELECTROLYSIS 

The results illustrated here have implications for catal- 
ysis, and electrolysis especially. In particular, we have 
shown that the hydrophobic surface that is formed from 
a passivated adlayer of water is accompanied by the 
existence of a liquid-vapor-like interface separating the 
adlayer from the bulk liquid. It has been previously 
demonstrated that excess protons [29] 130]. as well as some 
anions [31], preferentially adsorb to a liquid- vapor inter- 
face. It is expected, therefore, that the existence of this 
type of interface near the electrode surface itself acts cat- 
alytically by simply enhancing the local concentration of 
reactive species- protons - relative to the bulk. 

As one simple means of testing the assertion of a pro- 
ton enhancement at the interface, we have calculated 
the density profile for a fixed point charge model of 
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tion is calculated from, 

PH3O+O) = Pref e _/3F(z) 



(8) 



FIG. 7. Structure and density enhancement of hydronium 
ions at the liquid-vapor and electrode interface, (a) A char- 
acteristic snapshot of a HsO + at the liquid vapor interface, 
(b) A characteristic snapshot of a HsO + at the 111 inter- 
face, (c) Density distributions for the water (grey) and hydro- 
nium (blue), with a reference bulk density of the hydronium, 
Pref = 0.1g/cm 3 chosen for scale. 



a hydronium [32 cation using umbrella sampling. We 
choose this model as it has been shown previously that 
excess protons at the liquid vapor interface preferentially 
adopt a hydronium geometry over other forms, such as 
the Zundel cation [29] [30]. The local solvation structure 
of a hydronium in bulk is characterized by donation of a 
hydrogen bond by each of its three hydrogens, but its in- 
ability to accept any hydrogen bonds at the oxygen posi- 
tion due to the localized positive charge on the molecule. 
It has been demonstrated previously that this structure 
is conserved at the liquid vapor interface, and has been 
proposed as the reason for experimental observations of 
proton adsorption at the interface. A characteristic snap 
shot of such an expected configuration is shown in Fig. [7] 
(a). This snapshot is obtained from our simulations, in 
which find a mild enhancement of p(0) « 1.5pbuik a t the 
liquid- vapor interface, where the zero refers to the posi- 
tion of the Gibbs dividing surface. 

At the water-electrode interface the adlayer is com- 
posed almost entirely of hydrogen bond acceptor sites 
and thus it is expected that density of hydronium ions 
will be even further enhanced by their ability to donate 
hydrogen bonds into the adlayer. A characteristic snap 
shot of this type of configuration at the 111 surface is 
shown in Fig. [7](b). Figure [7] (c) confirms a much larger 
density enhancement of hydronium ions at the interface 
relative to the liquid- vapor case. This density distribu- 



where p re f is the density of the hydronium in the bulk 
and F(z) is the potential of mean force for moving the 
center of mass of the hydronium along the z direction. 
The enhancement found from this simplistic calculation 
is p(z*) ~ lOpbuik, much larger than the enhancement 
found at the liquid vapor interface. We note that while 
important aspects of proton derealization and polariz- 
ability are neglected in this calculation, each affect has 
been shown previously to further enhance interfacial ad- 
sorption. Thus while this calculation is meant only as 
illustrative of our findings, and not expected to be taken 
quantitatively, we nevertheless suspect this behavior to 
be conserved in more detailed models. 

In the specific case of hydrogen evolution at a platinum 
electrode, it is generally assumed that the reaction pro- 
ceeds through two steps: the Volmer step, where a proton 
is transferred form the bulk and discharged at the metal 
surface, followed by the Tafel step where two adsorbed 
protons combine to form hydrogen and desorb from the 
surface [33]. The latter is considered to be the rate deter- 
mining step. The enhancement of the proton concentra- 
tion at the interface is consistent with platinum's ability 
to easily transfer and accumulate protons on and near 
the surface, while the long time relaxation on the surface 
detailed here and elsewhere [20 undoubtably makes dif- 
fusion of adsorbents slow. While these results are consis- 
tent with mechanistic assertions for hydrogen evolution, 
further work must be done to explore the full implications 
of the effects illustrated here on catalysis. 



V. METHODS 

The system simulated consists of a slab of water in 
contact with an electrode on one side and with a free in- 
terface on the other side with a vacuum layer of 40 A. 
The electrode consists of three layers of atoms, totaling 
nearly 500 particles, held fixed in an FCC lattice with 
spacing of 3.92 A and with either the 100 or 111 facet 
exposed to the solution. A slab of water nearly 40 A 
thick was placed in contact with the electrode, and the 
dynamics of the nearly 1800 molecules are propagated 
using a Nose Hover integrator [34]], with SHAKE im- 
posing bond and angle constraints for the water as im- 
plemented in LAMMPS [35]. All simulations were run 
at 298 K. Periodic boundary conditions are employed in 
the plane parallel to the surface, and slab corrections to 
ewald summations are applied. Interactions between the 
water molecules are computed from the SPC/E potential 
[16] • The water electrode potential is modeled following 
Siepmann and Sprik [6] where the platinum water inter- 
action is a sum of two and three body terms, parame- 
terized to get the correct value of the adsorption energy 
and ground state geometry as determined by quantum 
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chemical calculations. Additionally, to model the polar- 
izable metal surface each electrode atom carries a Gaus- 
sian charge of fixed width but variable amplitude, which 
is updated at each timestep by minimizing the energy of 
the slab subject to a constraint of equal potential across 
the conductor. A more thorough description of the elec- 
trode model can be found in [6j [10] . 

The long relaxation times in the system relative to the 
timescales accessible by our simulations make ensuring 
equilibration on the surface difficult. In order to check 
the sensitivity of our results to their initial conditions 
we have prepared an ensemble of 40 independent sur- 
faces, for both the 100 and 111 crystal faces. Of the 
surfaces, 20 were produced through quenching the sys- 
tem at a rate of 10 K/ns from a system equilibrated 
at T = 400, and 20 from a process of vapor deposition 
where water molecules are exposed to the surface at a 
rate of 2 molecules/ 1 ns nm 2 , which corresponds closely 
to the equilibrium desorption rate. We have redone all 
of the calculations presented in the main text over this 
extended surface ensemble and found results that were 
indistinguishable from those presented above. 

The calculations of the excess hydronium ion were ac- 
complish using umbrella sampling along the z coordinate. 
For charge neutrality a small anion was placed in the wa- 
ter slab but kept at distances greater than 15 A from the 
hydronium. 



Appendix A: Surface Relaxation 

Relaxation on the electrode surface is slow, owing to 
the large ordered domains patterned by the underlying 
lattice symmetry. In order to quantify the timescales 
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FIG. 8. The correlation function for dipole orientations for 
waters on the 100 and 111 surfaces. 

associated with surface reorganization we calculate a wa- 
ter's dipole time correlation function given it starts and 
ends adsorbed to the surface, 

c MlM (t) = <&(*) • ft(o)Mi(t)]M*(o)]> , (Ai) 

where fa is the dipole vector for the i water molecule and 
hi is an indicator function which is equal to 1 if the center 
of mass of water molecule i is within 3 A of a electrode 
atom and is otherwise. 

Figure [8] reports on this function for both electrode 
geometries. Both have a characteristic time constants 
that decay over timescales longer than 1 ns. The long 
timescales associated with this relaxation make ensuring 
equilibration of the simulations difficult. In order to ad- 
dress this we have confirmed that all of our results are 
independent of the preparation of the initial condition. 
See methods for details. 
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